gisR

GIS-Review

Thomas Knecht, Ueli Mauch

11. Juli 2024

Was ist R

R ist eine offene Programmiersprache welche ursprünglich für statistische Zwecke erstellt wurde.

Der Einsatzbereich von R hat sich heutzutage weit über die statistischen Anwendungen etabliert.

Hinweis

Gerade für Data-Wrangling und Data-Engineering kommt die Stärke von R besonders zu tragen.

GIS mit R

Seit einiger Zeit hat sich R im Bereich GIS weiterentwickelt. Es lässt sich nun sehr einfach GIS-Prozesse mit Datenanalyse verbinden und es ist kein Medienbruch mehr nötig.

Packages (gleich wie libraries in python)

  • sf: Arbeiten mit Feature Classes

  • terra: Arbeiten mit Raster-Daten

  • stars: Arbeiten mit Spatio-Temporal Arrays

Einlesen der Daten

# Lesen der Gebäudedaten
# gwr <- read.csv("https://www.web.statistik.zh.ch/ogd/daten/ressourcen/KTZH_00002022_00004064.csv") |>
#   dplyr::rename_all(tolower)

gwr <- readRDS("gwr.RDS") |> 
  dplyr::rename_all(tolower)
#
# saveRDS(gwr, file = "gwr.RDS")

# Umwandeln in ein sf-Objekt -> Geodatenobjekt
gwr_sf <- sf::st_as_sf(gwr, coords = c("e.gebaeudekoordinate", "n.gebaeudekoordinate"), crs = 2056)

# interaktives Darstellen
mapview::mapview(head(gwr_sf, 100))

gemeinden <- sf::st_read("https://maps.zh.ch/wfs/OGDZHWFS?SERVICE=WFS&VERSION=2.0.0&Request=getfeature&TYPENAME=ms:ogd-0095_arv_basis_up_gemeinden_seen_f&outputformat=geojson") |>
  dplyr::rename_all(tolower)

Analysieren der Daten

# Matchen der Gemeindeinformation an die GWR-Daten
# Die geometrie wird hier entfernt, da die Datenanalyse danach viel schneller läuft
gwr_sf_gem <- gwr_sf |>
  dplyr::select(gebaeudekategorie_code, gebaeudekategorie_bezeichnung) |>
  sf::st_join(gemeinden) |>
  sf::st_drop_geometry()


# Berechnung der Anteile der Geböude ohne Wohnnutzung am gesamten Gebäudebestand einer Gemeinde
geb_ohne_wohnnutzung_pro_gem <- gwr_sf_gem |>
  dplyr::group_by(gebaeudekategorie_code, gebaeudekategorie_bezeichnung, bfs, gemeindename) |>
  dplyr::summarise(anzahl = dplyr::n()) |>
  dplyr::ungroup() |>
  dplyr::group_by(bfs, gemeindename) |>
  dplyr::mutate(anteil = round(anzahl/sum(anzahl)*100, 2)) |>
  dplyr::ungroup() |>
  dplyr::filter(gebaeudekategorie_code == 1060)

Vorbereiten für das Darstellen

# Hinzufügen der Gemeindepolygone
# Herausfiltern der Seen und fixen der Polygone
geb_ohne_wohnutzung_sf <- gemeinden  |>
  dplyr::select(bfs) |>
  dplyr::left_join(geb_ohne_wohnnutzung_pro_gem, by = "bfs") |>
  dplyr::filter(bfs != 0) |>
  sf::st_make_valid() |>
  # tmap vewendet die erste Spalte als Bezeichner, deshalb hier die Änderung in der Spaltenanordnung
  dplyr::select(gemeindename, dplyr::everything())


geb_ohne_wohnnutzung_sf_generalized <- geb_ohne_wohnutzung_sf |>
  rmapshaper::ms_simplify(keep = 0.005)

Statisches Darstellen der Daten

## Nur für Präsentation!
geb_ohne_wohnnutzung_sf_generalized <- readRDS(here::here("data.RDS"))

# Statisches Darstellung
tmap::tm_shape(geb_ohne_wohnnutzung_sf_generalized) +
    tmap::tm_polygons("anteil")

Interaktives Darstellen der Daten

## Nur für Präsentation!
geb_ohne_wohnnutzung_sf_generalized <- readRDS(here::here("data.RDS"))

# Setzt tmap auf interaktiv
tmap::tmap_mode("view")

# Interaktive Darstellung
tmap::tm_shape(geb_ohne_wohnnutzung_sf_generalized) +
    tmap::tm_polygons("anteil")

Remote Sensing mit R

Extensiv genutzte Wiesen müssen mindestens einmal pro Jahr gemäht werden und das Schnittgut muss abgeführt werden. Die Flächen dürfen in Abhängigkeit der Zone jeweils frühestens Mitte Juni bis Mitte Juli genutzt werden1

Fragestellung

Gibt es Landwirtschaftsflächen mit Schnittzeitpunkt 15. Juni, welche bereits vorher gemäht wurden?

Workflow:1

  • Sentinel-2 Daten herunterladen
  • NDVI berechnen
  • Statistik für Landwirtschaftsflächen errechnen

STAC

  • SpatioTemporal Asset Catalog1

  • Wir benutzen den Planetary Computer Data Catalog von Microsoft2

STAC

  • Collections abfragen
stac_source <- rstac::stac(
  "https://planetarycomputer.microsoft.com/api/stac/v1",
  force_version = "1.0.0"
)
collections_query <- stac_source |>
  rstac::collections()

rstac::get_request(collections_query)

STAC

# BBOX definieren
ktzh_bbox <- "8.3576933717031, \
              47.159435621727, \
              8.9849509769386, \
              47.694469737114"

# STAC Query bauen
stac_query <- rstac::stac_search(
  q = stac_source,
  collections = "sentinel-2-l2a",
  datetime = "2023-06-10/2023-06-14",
  bbox = ktzh_bbox
) |>
  rstac::ext_filter(
    `eo:cloud_cover` <= 20)

Filtern nach Gebiet, Zeitraum, Collection Sentinel-2-L2A. Bewölkung einschränken.

Mehr zum Download in scripts/dlFromSTAC.R

Geodaten von WFS laden

wfs_zh <- "https://maps.zh.ch/wfs/OGDZHWFS"

url <- httr::parse_url(wfs_zh)
url$query <- list(service = "WFS",
                  version = "2.0.0",
                  request = "GetFeature",
                  typename = "ms:ogd-0095_arv_basis_up_bezirke_f",
                  outputformat = "geojson")

request <- httr::build_url(url)
bezirke <- sf::read_sf(request)

# Save Bezirke Layer
sf::write_sf(bezirke, "geodata/Bezirke.gpkg")

WFS GetFeature Request, am Beispiel der OGD Bezirksgrenzen

Raster

  • Bänder-Kombination & Visualisierung (terra Package)
RGB FCIR
RGB FCIR

NDVI

Quelle

Rasterrechner

  • 4-Band Raster als Basis
  • NDVI nach dieser Formel

\[ NDVI = \frac{NIR - R}{NIR + R} \]

  • Berechnung simpel:
# calculate NDVI using the red (band 1) and nir (band 4) bands
sent2_ndvi <- (RGBI_zh[[4]] - RGBI_zh[[1]]) / (RGBI_zh[[4]] + RGBI_zh[[1]])

sent2_ndvi |> terra::writeRaster("geodata/ndvi_terra.tif")

Statistik

  • Nutzflächen auswählen, Feld area berechnen
  • 2 Filter anwenden
# Load Polygon Layer
LWNutz_Diels <- sf::read_sf("geodata/LWNutz_Dielsdorf.gpkg") |>
  mutate(area = sf::st_area(geom))

# Set aller Flächen, die erst am 15.06. sollten geschnitten werden:
BF_Fl <- LWNutz_Diels |>
  filter(harvest_date == "15.06.") |>
  # Drop die kleinsten Flächen unter 500 m^2
  filter(area > units::set_units(500, "m^2"))

Zonale Statistik

  • Statistik pro Nutzfläche, jeweils der mean wird berechnet
bfs_parz area mean_ndvi geom
0060_692 1534.220 [m^2] 0.4135512 POLYGON ((2681108 1262025, ...
0088_1041 2917.503 [m^2] 0.5237937 POLYGON ((2678319 1262664, ...
0088_965 10950.015 [m^2] 0.5349045 POLYGON ((2678276 1262340, ...
0058_3560 10104.811 [m^2] 0.5911278 POLYGON ((2678993 1268139, ...

Verteilung

Map

# Setzt tmap auf interaktiv
border <- c(2669379, 1252226, 2684597, 1269716)
breaks <- c(-0.4, -0.2, 0, 0.2, 0.4, 1)

tmap::tm_shape(BF_Fl, bbox = border) +
tmap::tm_polygons("mean_ndvi", style = "fixed", breaks = breaks) +
tmap::tm_shape(BF_Fl |> filter(mean_ndvi <= 0.4), bbox = border, name = "NDVI <= 0.4 ") +
tmap::tm_polygons("mean_ndvi", style = "fixed", breaks = breaks, legend.show = FALSE)
# TODO
tmap::tmap_mode("plot")
tmap::tmap_options(check.and.fix = TRUE)
border <- c(2669379, 1252226, 2684597, 1269716)
left <- tm_shape(BF_Fl,
                  bbox = border) +
tm_polygons("mean_ndvi",
            style = "fixed",
            breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 1)
            ) +
  tm_layout(title="ii")
right <- tm_shape(BF_Fl |> 
                    filter(mean_ndvi <= 0.4),
                  bbox = border) +
  tm_polygons("mean_ndvi",
              style = "fixed",
              breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 1)
  )
# TODO table
BF_Fl |> dplyr::filter(`mean_ndvi` <= 0.4) |> 
  dplyr::select(c("blw_name","region","mean_ndvi", "area")) |>
  sf::st_drop_geometry() |> head(10) |> 
  knitr::kable() |> kable_styling(font_size = 22)
blw_name region mean_ndvi area
Extensiv genutzte Wiesen (ohne Weiden) Chrüzstrasse 0.3896368 578.9319 [m^2]
Extensiv genutzte Wiesen (ohne Weiden) Chrüzacher 0.3629235 1498.1367 [m^2]
Extensiv genutzte Wiesen (ohne Weiden) Untergrüt 0.3985794 1536.4542 [m^2]
Extensiv genutzte Wiesen (ohne Weiden) Schöckfeld 0.3787448 2046.3655 [m^2]
Wenig intensiv genutzte Wiesen (ohne Weiden) Hard Humusdepot 0.3058146 10962.5786 [m^2]
Hecken, Feld-, Ufergehölze mit Krautsaum Riedthof Ost 0.3903728 798.4844 [m^2]
Extensiv genutzte Wiesen (ohne Weiden) Feisterloo 0.3868750 941.6327 [m^2]
Extensiv genutzte Wiesen (ohne Weiden) Lierehölzli 0.3605108 681.2977 [m^2]
Extensiv genutzte Wiesen (ohne Weiden) Untergrüt 0.3524443 4635.9250 [m^2]
Extensiv genutzte Wiesen (ohne Weiden) Fallmet 0.3872881 2092.0953 [m^2]

Fragen?